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We present Cahn-Hilliard and Allen-Cahn numerical integration algorithms that are uncondition- 
ally stable and so provide significantly faster accuracy-controlled simulation. Our stability analysis 
is based on Eyre's theorem and unconditional von Neumann stability analysis, both of which we 
present. Numerical tests confirm the accuracy of the von Neumann approach, which is straight- 
forward and should be widely applicable in phase-field modeling. We show that accuracy can be 
controlled with an unbounded time step At that grows with time t as At ~ t a . We develop a clas- 
sification scheme for the step exponent a and demonstrate that a class of simple linear algorithms 
gives a = 1/3. For this class the speed up relative to a fixed time step grows with the linear size of 
the system as N/ log N, and we estimate conservatively that an 8192 2 lattice can be integrated 300 
times faster than with the Euler method. 

PACS numbers: 64.75.+g, 05.10.-a, 02.60.Cb 



I. INTRODUCTION 

A starting point in the analysis of coarsening systems, 
such as the phase separation dynamics following a quench 
from a disordered to an ordered phase, is the characteri- 
zation of the asymptotic late-time behavior. Most coars- 
ening systems exhibit asymptotic dynamical scaling with 
the characteristic length scale L(t) given by the size of 
individual ordered domains. The growth-law L ~ f" is 
determined by only a few general features, such as con- 
servation laws and the nature of the order parameter (see 
PJ for a review). For conserved Cahn-Hilliard equations 
describing phase-separation, L ~ i 1 / 3 at late times. More 
detailed information about the scaling state is difficult to 
obtain analytically. Indeed the very existence of scaling 
has only been demonstrated empirically in simulations 
and experiments. Consequently, computer simulations 
of coarsening models, especially phase-field type mod- 
els like the Cahn-Hilliard equation, play an essential role 
in our understanding and characterization of late-stage 
coarsening. 

These simulations face several restrictions. To accu- 
rately resolve the asymptotic structure it is necessary to 
evolve until late times so that L(t) ^> w, where w is the 
domain wall width. However, to avoid finite-size effects 
we must halt the simulation when L(t) is some fraction 
of the system size L sys . Additionally, to resolve the do- 
main wall adequately the lattice spacing Ax must be be 
sufficiently small compared to the domain wall width w. 
Very large lattices of linear size L sys / Ax are necessary to 
satisfy all of these requirements: Ax < w <C L(t) < L sys . 
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Accurate studies of the scaling state require us to evolve 
large systems to late times. 

Unfortunately, current computational algorithms are 
very inefficient in their time integration. The standard 
Euler integration of the Cahn-Hilliard (CH) and Allen- 
Cahn (AC) coarsening models, for conserved and non- 
conserved dynamics, respectively, is known to be unsta- 
ble for time steps At above a threshold fixed by the 
lattice spacing Ax — this is the "checkerboard" insta- 
bility |2|. This imposes a fixed time step irrespective 
of the natural time scale set by the physical dynamics. 
The domain walls move increasingly slowly, for example, 
the CH equation yields asymptotic domain wall velocities 
v ~ dL/dt ~ t~ 2 / 3 . Consequently, a fixed time step re- 
sults in ever-decreasing amounts of domain wall motion 
per step and eventually becomes wastefully accurate. 

Ideally, one would like a stable integration algorithm, 
which would allow accuracy requirements rather than 
stability limitations to determine the integration step 
size. Recently, Eyre proved the existence of uncondition- 
ally gradient stable algorithms (essentially a strict non- 
increase in free energy for every possible time step) 0, 
and provided explicit examples of stable steps for both 
CH and AC dynamics 0, • The present work is con- 
cerned with developing these methods in two directions: 
clarifying and expanding the class of unconditionally sta- 
ble algorithms, and deriving the accuracy limitations on 
these algorithms. 

Our main results for stability are the following. 
We have determined the parameter range for which 
Eyre's theorem proves unconditional gradient stability 
(Sec. lII A)l . and we present Eyre's theorem in appcndixlAl 
We have also determined the parameter range that is un- 
conditionally von Neumann (vN) stable, that is, linearly 
stable for any size time step (Sec. Ill B]l . The latter range 
is a superset of the former, and neither appear to have 
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been previously determined. We have also performed nu- 
merical tests of stability in dimension d = 2 (Sec. Ill C[) 
and found that the vN stability condition appears to be 
sufficient for identifying unconditionally gradient stable 
steps. Specifically, for the parameterless form of the CH 
equation (see pj) 



-V 2 (V 2 + </>-0 3 ), 



(1) 



there exists a class of semi-implicit steps 

4> t+At + AtV 2 [(l-a 1 )4> t+At + (l-a 2 )V 2 4>t+At} 
= (f>t + AW 2 [-ai4> t -a 2 V 2 4> t + ^}. (2) 

that may be solved for the updated field <f>t+At efficiently 
by means of fast Fourier transform (FFT). The various 
stability conditions for these steps are depicted in terms 
of ai and a-i in Fig. ^ The stability conditions do not 
depend on the lattice type or dimension, on the volume 
fraction, or on the form of the lattice Laplacian. This 
implies, for example, that these algorithms could be com- 
bined with adaptive mesh techniques (see, for example, 
H) for independent control of spatial and temporal dis- 
cretization. Fig. suggests that the unconditional vN 
stability conditions, which are widely applicable and rel- 
atively easy to analyze, may provide a reasonably accu- 
rate proxy for unconditional gradient stability. We have 
also determined the analogous stability conditions for the 
AC equation. 
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FIG. 1: For time steps parametrized as in the dark shaded 
region indicates parameters for which Eyre's theorem proves 
unconditional gradient stability, while the light shaded region 
corresponds to unconditional von Neumann (linearly) stable 
steps. The open circles denote steps that are numerically 
gradient stable under all of our tests, as described in Sec. una 
while the black circles indicate parameters that were found 
numerically not to be gradient stable. 

When stability is not the limiting factor, practical lim- 
its are still imposed by accuracy. To maintain the domain 
wall profile to a given accuracy, a time step should be 
chosen so that the wall only moves a fraction of its width 
w in a single step. For a scaling system with L ~ t n , 
where n < 1 generally, the passage time r scales like 



w/v ~ w/L ~ t 1 "at late times, 
step should scale as 

Ai„„f ~ t ~ t 1 



Then the natural time 



(3) 



For CH dynamics, n = 1/3 and At nat ~ t 2 / 3 while for 
AC dynamics n = 1/2 and At na t ~ t 1 / 2 . However, we 
show that these stable algorithms arc still not capable of 
accurately simulating coarsening using the natural time 
scale — despite their stability. For example, accuracy 
limits the stable CH steps given above to "only" At ~ 

To understand the limitations imposed on even stable 
algorithms by accuracy, we study in Sec. 11111 the trunca- 
tion error for the CH equation for general numerical al- 
gorithms, and determine the how these terms scale with 
time to all orders in At (Sec. IIII D|l . We develop a classi- 
fication scheme for such algorithms based on the lowest 
order p of At p at which truncation error fails to follow 
its optimal scaling and show that this term limits the ac- 
curacy of the algorithm at late times (Sec. IIII A|) . Our 
analysis leads to the conclusion that accuracy requires a 
time step 



At „ t 2( P -l)/3p 



(4) 



for the CH model. The algorithms in Eq. (j5J have p = 2, 
meaning the error becomes sub-optimal at 0(Ai 2 ), the 
leading error term. This result is consistent with our 
numerical observations. Our simple analysis for the nat- 
ural time step, Eq. (J2J), corresponds to the p — oo class. 
We are unable to identify any such "perfect" algorithms 
for the CH case; they are quite likely impossible for any 
nonlinear problem. 

Next, we turn to the question of practical advantage. 
Various computational algorithms have been developed 
to mitigate the impact of instabilities by increasing At by 
a fixed factor compared to the simplest Euler discretiza- 
tion. For example, the cell-dynamical-scheme (CDS) 
exploits universality to choose a free energy that is con- 
venient in terms of numerical stability. More recently, 
Fourier spectral methods 0, have been shown to in- 
crease the maximum At by an impressive two orders of 
magnitude. However, these methods still require fixed 
time steps and so cannot adjust to the naturally slowing 
CH dynamics. 

In Sec. II VI we determine the relative advantage of in- 
tegration by algorithms such as Eq. compared to the 
conventional Euler method. For a reasonably conser- 
vative choice of accuracy requirements, we find for an 
8192 x 8192 lattice (currently feasible for a linux work- 
station) with Ax = 1 that the new methods can integrate 
up to finite size effects roughly a factor of 300 times faster 
than possible with the Euler method. The advantage of 
unconditionally stable steps increases with larger system 
sizes: for lattices of linear size N we show the relative 
advantage in speed is order N/\ogN, regardless of spa- 
tial dimension of the system. This means that as com- 
putational power continues to increase, unconditionally 
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gradient stable algorithms will become even more valu- 
able. 

We present a summary and outlook for future devel- 
opments and applications in Sec. [V] 



II. STABILITY 

The parameterless form of the CH equation for a con- 
served scalar field Q is 

= V 2 M (5) 
where /i is the local chemical potential given by 



and F[<fi] is the free energy functional, taken here to be 

F[4>]= J ^[i(V0) 2 + i(0 2 -l) 2 ]- (?) 

The second term in F represents a double-well potential 
with equilibrium values (fa = ±1, and Eqs. ©, and 
combine to give Eq. QJ. The parameterless form of 
the AC equation pj is 

0= — A*= V 2 cf> + 4>-cl) 3 . (8) 

For dissipative dynamics such as the CH and AC equa- 
tions, a discrete time stepping algorithm is defined to be 
gradient stable only if the free energy is non-increasing, 
F[<t>t+At] < -F[0t]; f° r an y field configuration <p t - The 
other requirements for gradient stability, e.g. that sta- 
ble fixed points must correspond to minima of F, or 
that F should increase without bound for large (j), are 
already manifest in the discretized forms of these equa- 
tions. Gradient stability may reasonably be regarded as 
the ultimate stability criterion for the CH equation. 

Unconditional gradient stability means that the 
conditions for gradient stability hold for any size time 
step At £ [0, oo). Since unconditionally stable steps 
are our primary concern, we will henceforth use "sta- 
ble" or "unstable" to refer to the behavior for arbitrarily 
large At. That is, "stable" implies unconditionally sta- 
ble, while a fixed time step algorithm like the Euler step 
may be referred to as "unstable" or conditionally stable. 

The Euler time discretization of the CH equation is 

0f + u At = t + AiVV (9) 

The Euler update is "explicit" since the field at the earlier 
time step {<pt) explicitly determines the field at the next 
time step (4>t+At)- It is also unstable for values of At that 
exceed a lattice-dependent threshold, At max ~ Ax 4 0- 
The fully implicit time step is obtained by replacing /i t 
with fit+At m Eq. ©, and is, like the Euler step, accurate 
to O(At). Other time steps, which involve splitting [i into 



parts evaluated at t and at t + At, are generally called 
semi-implicit methods. 

Remarkably, Eyre 0, 3 proved that appropriate semi- 
implicit parametrizations can lead to stable update steps 
for both the CH and AC equations. To explore these 
possibilities, it is useful to introduce a general family of 
such steps for the CH equation in an arbitrary spatial 
dimension: 

&+At + (l-a 1 )AtV 2 j> t+At + (l-a 2 )AtV l 4> t+At 

(l-a 3 )AtV 2 [Cfe"] = 
<j>t - aiAtV 2 t - a 2 AtV% +a 3 AtV 2 ^. (10) 

This reduces to Eq. J2J for a 3 = 1 . For each of the three 
terms on the right-hand side of Eq. there generally 
are both explicit and implicit contributions to Eq. llOj l. 
and this will be exploited to construct stable dynamics 
for any size At. For all values of the parameters and m 
this step gives a solution 4>t+At that is order O(At) accu- 
rate. The implicit terms are denoted 4>t+At, with cftt+At 
reserved to represent the exact field obtained by integra- 
tion of Eq. over the time step At. We choose our pa- 
rameterization such that a\ = a 2 = a 3 = 1 corresponds 
to the Euler update Eq. ©, while a\ = a 2 — a 3 = is 
the fully implicit step. For a 3 =/= 1 we have, motivated 
by Eyre, a mixed non-linear term with < m < 3 that 
combines implicit and explicit terms. 

It is useful to sort algorithms described by Eq. I|1(J|) 
into three categories based on how they are implemented 
numerically. First, when a 3 = 1 we have linear di- 
rect steps, where the equation for (f>t+At is linear and 
has spatially uniform coefficients so the updated field 
can be found efficiently with FFT methods. Second, 
when a 3 7^ 1 but m — 2 then the implicit equation re- 
mains linear in 4>t+At but no longer has spatially uni- 
form coefficients. Eyre outlines an iterative procedure 
for solving these equations [if], so we call these linear 
iterative steps. Insisting on convergence of the itera- 
tive procedure restricts this class to a subset of param- 
eter values. Finally, for a 3 =/= 1 and m =/= 2 the up- 
date equation is nonlinear. For some parameter values 
the nonlinear equation can lead, unphysically, to multi- 
ple solutions. This occurs for both the fully implicit case 
ai = a 2 = a 3 = m = 0, as well as the Crank-Nicholson 
case a\ — a 2 = a 3 = 1/2, m — 0, whenever At exceeds a 
threshold value Q . Generally the nonlinear equations re- 
quire solution by the Newton-Raphson method, which is 
complicated to implement in two or more spatial dimen- 
sions. For some parameter values this can be demon- 
strated to be absolutely convergent, so nonlinear steps 
provide a viable option — though not one we have ex- 
plored numerically. 

The step parametrization for the AC equation analo- 
gous to Eq. (fTUft is 

4>t+At - (1 - ai)At0 t+At - (1 - a 2 )AtV 2 4> t +At 

+ (l-a 3 )At[Cfe] = 
4> t + a-iAtfo + a 2 AtV 2 (j) t - a 3 At$, (11) 
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which we include because the theoretical stability anal- 
ysis follows nearly identically for the CH and AC equa- 
tions, and the stability regions are given by the same 
shaded regions of Fig. 1. 

A. Unconditionally Stable Steps from Eyre's 
Theorem 

Eyre's theorem (see appendix EJ shows that an un- 
conditionally gradient stable algorithm results, for both 
the CH and AC equations, if one can split the free en- 
ergy appropriately into contractive and expansive parts, 
F = F c + F E , and treat the contractive parts implic- 
itly and the expansive parts explicitly. That is, the CH 
equation J5J is discretized as 

4> t+At - AtV^t+At =4H + AiVVf , (12) 
while the AC equation © is discretized as 

4> t+A t + A^t+At =4t- At/if, (13) 

where fif — dF x /d<j>i for lattice site i, and where V 2 
implies a lattice laplacian. The necessary condition on 
the splitting is the same for both equations and may be 
stated by introducing the Hessian matrices 

d2F *,E 92 PE »rC 9 2 F C 

M- ■ = M = M — 

(14) 

where i,j denote lattice sites. First, we must have all 
eigenvalues of M. E non-positive and all eigenvalues of 
M non-negative. Second, as shown in appendix [S] for 
A m m equal to the smallest eigenvalue of M and A^ aa , the 
largest eigenvalue of M £ , we need 

^max — 2 ^min- 0-5) 

This also automatically satisfies the convexity require- 
ment for 'Ml E , since X m in < 0. 

To identify the appropriate splittings, it is useful to 
break the free energy Eq. JJJ, in its lattice-discretized 
form, into three parts (neglecting the irrelevant constant 
V/4 term): 

i i 

^ (3) = EM ■ ( 16 ) 

i 

with corresponding Hessian matrices M^K The first, 
= —Sij, where §ij is the Kronecker ^-function, has 

all eigenvalues equal to — 1. Next, M-^ — (— V 2 )^ is 
negative the lattice laplacian, which can always be diago- 
nalized by going to Fourier space. It immediately follows 
that the eigenvalues of M^ 2 ) are strictly non-negative. 
(Even for irregular spatial discretizations, the M^ 2 ) eigen- 
values must be non-ncg ative.) Finally, M^f = 3^6^, 



which has strictly non-negative eigenvalues as well. We 
parametrize the splitting via 

3 3 

F E = J2 a i Fii) F c = Y,( 1 ~a t )F^ (17) 

»=i i=i 

which results in the general CH step Eq. and AC 
step Eq. when m = 0. 

Now to obtain bounds: since the sum of matrices, M = 
M^ 1 ) +M( 2 ) +M( 3 \ has eigenvalues bounded by the sum 
of the bounds, the minimum eigenvalue of M satisfies 
Amm > — 1. Therefore Eq. I|15|) is satisfied by ensuring 
X E < -1/2. 

One example that satisfies these conditions is the split- 
ting F E = and F c = F^ + F<- 3 \ since \ E ax = -1 
satisfies Eq. (II 51) and M c has strictly non-negative eigen- 
values. This provides a gradient stable nonlinear step 
with a± = 1 and a 2 =03 = 0. This case was identified 
by Eyre |3(, who noted that the convexity requirement 
for the splitting guarantees absolute convergence of the 
Newton-Raphson method. 

Eyre also presents a technique for identifying stable 
linear direct algorithms Q , which relies on the fact that 
cj) 2 is bounded. It exceeds unity only slightly in the CH 
equation and only in the interior region of a curved in- 
terface due to Gibbs-Thompson effects 9] . Therefore the 
eigenvalues of M' 3 ) have an effective upper bound, ap- 
proximately three. If we then take F E — aiF^ + F^ 
(so as — I and a 2 = 0) the eigenvalues of M B are of the 
form -ai+3(/) 2 and satisfy Eq. (JTSJl for <j)f < 1 if a x > 7/2. 
Any value a 2 < will give the same result, since nega- 
tive values of a2 can only decrease the eigenvalues of M £ . 
These choices imply F c = (1 - a^F^ + (1 - a 2 )F^ 2 \ 
which has the necessary non-negative eigenvalues for the 
range of a\ and a 2 given above. Therefore we can iden- 
tify a class of gradient stable direct CH and AC steps 
as 

01 > 7/2 a 2 < a 3 = 1. (18) 

This gives the dark gray shaded region in Fig. ^ These 
represent sufficient restrictions on the to satisfy the 
conditions for Eyre's theorem; however other values of 
the ai may be gradient stable as well. 

Eyre provided specific step examples for all three im- 
plementation categories: a nonlinear step, with a\ = 1, 
o-i = a 3 = 0, and m = 0, a linear iterative step with 
772 = 2 and the same as the nonlinear step, and a lin- 
ear direct step with a\ = 3, a 2 = 0, and CI3 = 1 0, Q. 
The nonlinear step is the example presented earlier in 
this section, and its gradient stability follows from Eyre's 
theorem. However, it is not clear to us that Eyre's the- 
orem can be directly applied to the iterative steps, and 
in fact we find Eyre's proposed iterative method to be 
numerically unstable, as described in Sec. Ill CI Finally, 
the ai value in the direct step violates Eq. Q1H|). so this 
case does not follow from Eyre's theorem. 
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B. Unconditional von Neumann Stability 

Since Eyre's theorem provides, in principle, only a sub- 
set of the possible gradient stable steps, complementary 
approaches for determining stability are desirable. In 
this section we extend von Neumann's (vN) linear sta- 
bility analysis 10] to arbitrary time steps, which we call 
unconditional vN stability. Since any gradient stable al- 
gorithm is likely also linearly stable, the von Neumann 
analysis would appear to identify a superset of possibly 
gradient stable algorithms: in principle the vN analysis 
could also identify some unwanted non-linearly unstable 
algorithms. As shown in Fig. ^ though, the vN stability 
boundary corresponds quite well with the numerically de- 
termined gradient stability line. This leads us to suggest 
that the approach of imposing unconditional vN stabil- 
ity on a broadly parametrized class of semi-implicit al- 
gorithms, followed by numerical checking, could be fruit- 
fully adapted to a wide variety of applications. 

We analyze the general step Eq. (|1(J|> for linear sta- 
bility around a constant phase <p — <pQ. It is important 
to realize there are physical, and therefore desirable, lin- 
ear instabilities in the continuum CH and AC equations. 
Therefore it is important to distinguish between these 
and the unphysical instabilities induced by the numerical 
implementation. Take 0(x, t) = </>o+?7(x, i), and linearize 
the CH equation JTJ in 77 to get 77 = — V 2 (V 2 ?7+r7— 30q?7). 
Fourier transform this to get 

Vk = -A k (A k + k%)rik, (19) 
k 2 = 1-30 2 . (20) 

Here Ak is the eigenvalue of the Laplacian and is non- 
positive, with Xkmin < A k < (note that Ak = —k 2 
in the continuum). The minimum value Xkmin depends 
on the lattice, spatial dimension, and specific form of 
the laplacian. Similarly, for the same <fr linearize the AC 
equation |JS} in 77 and Fourier transform to get 

= (A k + fco)?7k. (21) 

The physical instability for both Eqs. i|19|) and H21|) oc- 
curs for 

- A k < kl (22) 

which corresponds in the CH equation to spinodal decom- 
position pj . We stress that while these Fourier modes are 
linearly unstable, the dynamics of spinodal decomposi- 
tion is gradient stable and represents a physical decrease 
of the free energy, which is why it must be retained. 

We now linearize and Fourier transform our general 
CH step Eq. I|10|) as above to get 

[1-A U A*{ (01 - 1) - A k (l - 02) 

+00 (1 - a 3)(3 - m)}]?7k,t+At = 

[1 - A k Ai{ ax + A k a 2 

+0oMa3 + m(a 3 - l))}]ry M (23) 



Writing this as 

[1 + AtC^t+At = [1 + At7£ k ]r7k, t , (24) 

the von Neumann stability criterion is 

|1 + AtC k \ > |! + Atftkl, (25) 

which implies that small deviations from the constant 
solution evolve to decrease in magnitude. We want to 
impose this stability condition for all k and arbitrary pos- 
itive At. For large At, Eq. (gSJ) implies |£ k | > \K k \. The 
left-hand side of Eq. (|2*5|) can be made to violate the in- 
equality for small At unless £ k > 0. Combining these 
conditions we have 

£k > \Hk\, (26) 

which is a necessary and sufficient condition for uncondi- 
tional linear stability. This condition applies to all first- 
order time steps that can be expressed in the form given 
by Eq. QZQ. 

We examine the linear stability condition in two steps. 
First, £ k >K k : 

0<£ k -^ k = (-Ak)[-l-A k + 3^]. (27) 

This reduces to the spinodal condition, Eq. Q22[l. Note 
that all the parameters (ax, 0,2,0,3,711) are absent from 
Eq. H27|) . so we cannot interfere with the spinodal con- 
dition. This evidently follows from having a first-order 
accurate step. Next, we check for £ k > — 7^ k , which gives 

2ai-l-[(3-m)(2a 3 -l)+m]0o + A k (2a 2 -l) > 0. (28) 

If we choose a 2 < 1/2, then since A k < we get 2eti — 1 — 
[(3 — m)(2a 3 — 1) + to]0q > 0. For a 2 > 1/2 we obtain a 
lattice-dependent condition, that is, our inequality would 
contain X kmin . 

We choose to restrict ourselves to lattice- independent 
stability conditions as these are more practical: they 
carry over into any lattice or spatial dimension. For this 
purpose we take a 2 < 1/2. This gives the vN stable 
conditions 

a 2 < 1/2 

1 + max[0, (3 - m)(2a 3 - 1) + to] 
a x > . (29) 

We have let 0q vary in the late-time asymptotic range of 
0o S [0, 1], where Gibbs-Thompson induced supersatura- 
tion has be ignored, and have imposed on ax the most 
restrictive value that results. For this reason algorithms 
near the stability boundaries should be avoided at early 
times. 

For direct steps, with 03 = 1, the second condition in 
Eq. (|29() becomes ax > 2. This gives the lightly shaded 
region in Fig. ^ The Euler update, with ai = a 2 = 
as = 1 is clearly unstable since a 2 > 1/2 and a\ < 2. For 
linear iterative steps, with to = 2, Eq. (|29|l becomes ax > 
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max[l/2,a3 + 1]. The stability condition of the general 
nonlinear step cannot be further simplified from Eq. (|29[l . 
but the special case m = gives a± > max[l/2,3ci3 — 1]. 

There is another special case for which the stability 
conditions can be imposed, namely when m = and 
a± = a,2 = as = a. In this case the vN condition Eq. I|28|l 
becomes 

(l-2o)[-l + 3^ + Ak] >0. (30) 

The square brackets term is again the spinodal condition 
and should be positive for all physically stable modes, so 
for a < 1/2 both vN stability conditions reduce to the 
spinodal condition. However, these steps, which include 
the marginal Crank-Nicholson case (a = 1/2) and the sta- 
ble fully implicit step (a = 0) suffer from having multiple 
solutions to the nonlinear implicit equation whenever At 
exceeds some threshold, making them unsuitable. 

Regarding Eyre's proposed steps, introduced at the 
end of Sec. Ill Al we note that the direct step is vN sta- 
ble, the iterative step is marginal for vN stability, and 
the nonlinear step, which was gradient stable by Eyre's 
theorem, is also vN stable. 

The same linearization for the general AC step Eq. i|ll|) 
results in the same linearized equation ill' Ml) but with the 
substitution — X^At — > At. Since — Ak > 0, the vN sta- 
bility analysis for the AC equation is identical and also 
results in Eq. 

C. Numerical Stability Tests 

The vN stability analysis yields a considerably larger 
parameter range for stable steps, Eq. (|29|l . than those 
which are provably stable by Eyre's theorem, e.g. 
Eq. i|18|) . Here we determine numerically which step 
parametrizations are gradient stable, for purposes of 
comparison with the theoretical results. We focus pri- 
marily on direct steps, with as = 1, since these are an 
important practical class of steps. We consider only sym- 
metric quenches of the CH equation in this section, with 

W = o. _ 

The primary result, shown in Fig.^ is obtained as fol- 
lows. We evolved a uniformly distributed 20 x 20 array of 
direct CH steps with the parameter values a\ € (0, 4) and 
a-2 € (—2,2) on a 512 2 lattice to a final time t max . We 
take lattice spacing Ax = 1 here and throughout. At reg- 
ular intervals during the evolution we tested a single di- 
rect step with < At < 10 10 . This step was only used for 
stability testing, and did not contribute to the time evo- 
lution. Steps larger than At = 10 10 were not employed, 
to avoid spurious roundoff error effects. Any system that 
ever increased its free energy was labeled unstable, and 
plotted in Fig. ^ with a filled circle. The systems were 
evolved in time with multiple methods. First, we used 
Euler updates (At = 0.05) evolved to t max = 10 4 . Next, 
we evolved systems with direct updates both with fixed 
At = 100 and with an increasing time step At = 0.05 i 1 / 3 
(both to t max = 10 6 ). 



As Fig. ^ shows, all vN stable algorithms were found 
numerically to be gradient stable, and the lightly shaded 
region corresponds extremely well to the gradient sta- 
ble systems. Indeed, the vN stability boundary for ai 
appears to be followed quite sharply in the numerical 
tests. We do find numerical gradient stability for a re- 
gion where a 2 > 1/2: this is most likely due, ironically, 
to a lattice-induced stabilization. That is, since the lat- 
tice laplacian Ak has an implementation-dependent mini- 
mum value, the inequality i|28f) may be satisfied for some 
a2 > 1/2. Therefore we expect the precise location of 
this boundary to shift slightly depending on the lattice, 
the spatial dimension, and the choice of lattice laplacian, 
but not to cross 02 = 1/2. 

With the numerical tests described above we have 
tested the linear iterative step proposed by Eyre Q and 
found it to be unstable. 

To help illustrate numerical testing of gradient stabil- 
ity, we show a mixture of stable and unstable steps in 
Figs.0and[21 The difference between gradient stable and 
unstable steps is striking: while neither are particularly 
accurate for extremely large At, the unstable steps show 
a marked increase in the free energy density, while the 
gradient stable steps adhere to the strict non-increasing 
free energy condition. However, the closer view in Fig. 
shows that some cost is paid in accuracy: for small values 
of At, both the Euler step and the unstable semi-implicit 
step track the physical behavior better than the stable 
step. While it may appear from Fig. |3| that moderately 
large steps may be used with unstable algorithms, this 
is not case: for example using a At > 0.05 for the Euler 
update will lead to instability via accumulated error from 
repeated steps. 



III. ACCURACY 

With a gradient stable algorithm, it is possible to use 
a progressively larger time step as the characteristic dy- 
namics become slower. The limiting factor for the in- 
crease of the time step is then an accuracy requirement. 

A specified accuracy criterion may be imposed on the 
stable steps identified in Sec. UTI without any further the- 
oretical development using standard numerical adaptive 
step-size techniques (as described in [lOj and discussed 
in Sec. IIII Bll . Naively, one would expect a time step 
growing as At ~ t 2 / 3 , for the reasons presented in Sec. [I] 
However, this is not the case: empirically we find signif- 
icantly slower growth. This motivated us to study the 
sources of error terms in the gradient stable steps. Our 
main result is the p classification scheme, which deter- 
mines the allowed growth rate of the time step according 
to Eq. 0. 
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FIG. 2: Plot of the free energy density e versus time (thick 
solid line) approaching the asymptotic e ~ t -1 / 3 decay, as 
evolved with a Euler update with At = 0.01 in a 1024 2 system. 
At five distinct departure times td, separated by factors of 4, 
we show the free energies that result from a single time step 
At e (0,10000), plotted versus t = t d + At. The dotted 
lines correspond to using a common semi-implicit algorithm 
(oi — 1, a2 — 0, a$ = 1) for the single step, while the thin 
solid lines correspond to single steps with a vN stable direct 
algorithm (ai = 3, a2 = 0, and as = 1). 
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FIG. 3: As per Fig. H but with t d = 1000. The dashed 
line corresponds to a single step of the Euler update, which 
is gradient unstable. Both the Euler step and the unstable 
semi-implicit step (dotted) are unstable under repeated steps 
for much smaller At than appear to be accurate for a single 
step. 



A. The p Classification Scheme 

We begin with an analysis of the error magnitude asso- 
ciated with the various gradient stable algorithms. The 
exact tfit+At, obtained by integration of Eq. from a 
given <fi t , can be expressed in terms of the fields at time 
t by means of a Taylor expansion: 



(fit+At 



\At 2 d 2 ^ 



±At 3 d 3 < 



(31) 



The Euler update, Eq. (|5J), is simply the truncation of 
this expansion at O(At) with resulting error A(fi Eu = 



D t+At 



k+At given by 



A(f> 



Eu 



At n 



71=2 



d?(fit. 



(32) 



Other step parametrizations will have different coeffi- 
cients for the 0(At n ) component of the error, but the 
general feature of an expansion to all powers of At will 
be the same. Since our goal is to have a growing time 
step with controlled error, successively higher powers of 
At will require coefficients decaying increasingly faster in 
time. In order to determine the limitation on how fast 
the time step may grow, it is essential to know the decay 
rates of the coefficients of At n to all orders n. In this 
section we demonstrate how this can be done. We make 
use of the following results for asymptotic decay rates, 
derived in Sec. lIIIDl In the interfacial region (defined in 

sec.inni 



t -2n/3 «9«(v 2 ) V' ~ i~ 2n/3 (33) 



whereas in the bulk, that is, all of the system not near 
an interface, we find 

dt(fi ~ ^-( 1 /3)-(2/3)« Qn(y2\k^ ^ t -(l/3)-2(n+fc)/3 

(34) 

Consider first the Euler step: all the 0(At") coeffi- 
cients are simply proportional to the time derivative d 7 t l (fi 
evaluated at t. If numerical stability were not a prob- 
lem and we simply increased the time step according to 
the naive At ~ i 2 / 3 , we would find in the interfacial re- 
gion that every order in the Taylor expansion provides an 
O(t ) contribution to the error, whereas in the bulk re- 
gion every order provides an 0(t -1 / 3 ) contribution. This 
would present an accurate solution with a At ~ i 2 / 3 time 
step, except that of course the Euler step is not gradient 
stable for large time steps. 

Now consider the general step, Eq. (|10|l . The error 
term in this step, A<fi = (fit+At — <fit+At can be written as 



A4> = A(fi Eu 

-(l-ai)AiV 2 (0 t+A t-^t) 
-{\-a 2 )AtV i ^ t +At-(fit) 
+ (l- a3 )AiV 2 [C(^A t - 



4~ m )\ (35) 



This peculiar form with implicit (fit+At on the right is 
useful for the error analysis. By using Eq. (|10|) iteratively, 
the implicit terms can be replaced by terms higher order 
in At involving the field <fi t . For example, we can derive 
the 0(At 2 ) part of the error, using (fit+At — (fit = At (fit + 
0(Ai 2 ) and 4>lXAt-$~ m = {3~m)At(l> 2 t - m <fit + 0(At 2 ). 
We find the error in our general step to be 



A<fi 



\<fit + (<2l 



1)V 2 ^ + (a 2 - l)V 4 <fi\ 



+ (l-a 3 )(3-m)V 2 



• J t H>t 



At 2 + 0(At 3 )(36) 
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where the first term comes from Eq. ().'>2t . Now compare 
the time decay of the various terms. At the interface, the 
4>t part decays as t~ 4 / 3 , but the other terms all decay as 
£-2/3 Therefore, for general values of the ai and to, to 
keep the 0(At 2 ) interfacial error fixed the time step is 
limited to stow as At ~ t 1 ' 3 . We sec that the Euler case 
was special because it made all but the first term in the 
0(At 2 ) error vanish. Since every term in Eq. (|36[l decays 
faster in the bulk than at the interface, we conclude the 
error is interface limited, i.e., the accuracy criterion at 
the interface will determine how fast the time step can 
grow. This is a generic feature, as we will show below. 

There are other ways besides using the Euler step to 
make the 0(At 2 ) interfacial error decay as t~ 4 / 3 . If the 
coefficients satisfy 

a\ = a-2 = 1 — b = 1 — 36/ (3 — to). (37) 

for some 6, then the various <j>t terms in Eq. (|36|) add to 
give b(f>t- In this case, 

A<f> = -Ai 2 ^-b)4> t At 2 + 0(At 3 ) (38) 

and so the 0{At 2 ) coefficient decays as i~ 4 / 3 at the in- 
terface, and faster in the bulk. From this example we can 
construct the p classification scheme. 

Consider the truncation error term of order At n . This 
can be obtained by iterating Eq. Ij35(l and can be ex- 
pressed as a sum of terms of the form 9™~ 1 (V 2 ) fc J . If 
these terms appear in the right proportions, they com- 
bined via Eq. m to become proportional to <9™</>, which 
decays faster by a factor of 1/t 2 / 3 at the interface. This 
is exactly what occurs in the n = 2 case above when 
Eq. J23 is satisfied. 

Now consider some value p > 2 for which all At n 
error terms with n < p are proportional to At n d"(f>t, 
but at order to > p this breaks down into a sum of 
terms of the type At m d^~ l (V 2 ) k c^ t . In this case the 
order p term provides the leading asymptotic error. Fo- 
cusing on interfacial region, the order p term goes as 
^P£-2(p-i)/3 accorc ii n g to the second term in Eq. (J33J). 
Choosing the time step to hold this term at constant er- 
ror would require At ~ t a with a = 2(p — l)/(3p), as 
displayed in Eq. @. Now we show that all higher- and 
lower-order terms in At will decay faster than the At p 
term for this choice of a. For n < p, we have from the 
first term in Eq. (33J At n t- 2n / 3 ~ f l ("-2/3) _ t -2n/z P ^ 
so the n < p terms give ever-decreasing contributions 
to the error. For to > p the error terms are of the 
form At m t- 2 ( m -^/ 3 ~ i-2(m-p)/(3p) which decay as welL 

Hence the asymptotic interfacial error is given by the 
0(At p ) term as advertised, and is order t . Note that 
for this interface limited At ~ t a all bulk terms to all 
orders have decaying error terms, thus establishing inter- 
face limited error as a generic feature. 



B. Quantifying Error for Direct Steps 

Direct steps, with 03 = 1 by definition and a\ > 2, 
a 2 < 1/2 for stability, fail to satisfy Eq. (|3TJl and so all 
direct steps give p = 2 algorithms with At ~ t 1 ^ 3 . This 
means that the asymptotic error magnitude should be 
given exactly by 

|A0| = At 2 \( ai - 1)V 2 + (02 - 1)V 4 0| (39) 

with At — At 1 ' 3 . This gives a fixed amount of error at 
the interface, and all higher orders of At give decaying 
contributions. Therefore, the error magnitude is propor- 
tional to A 2 , and we can use numerical measurements of 
Eq. H39|) to develop the constant of proportionality. 

We determine error numerically in the usual way |1C| : 
compare the field 4>^ obtained from a single step of size 
At to the field 4>^ obtained from two steps of size At/2. 
It is straightforward to show that if the true error of the 
step is EAt 2 + 0{At 3 ), then 0« - ^ = (E/2)At 2 + 
0(At 3 ). Since we expect exactly At 2 error, we simply 
take 2[(fP-> — 4>^) to be the true error. 

In the bulk, the error decays as t~ 2 ^ 3 . The inter- 
facial error is not decaying, but the amount of inter- 
face decays as t^ 1 / 3 , which means the error magnitude 
Eq. I|39|) averaged over the entire system will also de- 
cay as i -1 / 3 , all from the interfacial contribution. To 
determine the error per lattice site in the interfacial re- 
gion, it is necessary to divide the averaged error by the 
fraction of the system in the interfacial region. We do 
that as follows. The asymptotic free energy density is 
given by the product of the surface tension a and inter- 
face density: e(t) = a A int (t) / 'L d sys ~ t" 1 / 3 , where the 
interfacial "area" Ai nt is a d — 1 dimensional hypersur- 
face, and L sys is the system size. For interface width 
w, A int (t)w/Lj ys — we /a represents the fraction of the 
system in the interfacial region. Multiplying the aver- 
aged error by a /(we) then gives the typical error in the 
interfacial region. The surface tension corresponding to 
Eq. (|7|) is a = 2^2/3. We take w = 2^2 as a typical 
measure for the interface width. 

We have investigated this error for a variety of direct 
algorithms in Fig.0J where we have plotted the interfacial 
error as determined above divided by A 2 . We plot this 
error amplitude against a\ and ai for the same shaded 
regions ["vN" and "E"] as identified in Fig-QJ The typical 
interfacial error for a given direct step of size At = At 1 / 3 
may be obtained by multiplying the appropriate contour 
value by A 2 . 

To illustrate the advantages of stable algorithms, as 
well as of a detailed error analysis where it is possible, 
we show in Fig. [5] how the error evolves in time for di- 
rect steps with At = At 1 / 3 versus the Euler step with 
fixed At. The field <j> is evolved by the Euler method, 
and during the evolution error checking is done with sin- 
gle steps that do not contribute to the evolution. The 
decay of the Euler error shows that the Euler method is 
asymptotically wastefully accurate. 
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FIG. 4: Contour of scaled error for a single direct update in a 
1024 2 system. The systems are evolved well into the scaling 
regime it ~ 3000) with a fixed-step Euler update. The errors 
are found by comparing a single direct time step At = At 1 ^ 3 
with two steps of size At/2, and are then scaled by 2a/(A 2 we) 
to estimate the average error magnitude per lattice site in the 
interfacial region, as described in the text. 
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FIG. 6: Plot of e versus t for a Euler update (with At = 0.05, 
thick solid line) and with the evolution via a direct algorithm 
(ai = 3 and a 2 = 0) driven with At = At 1/3 with A = 0.1 
(dotted line) and 0.01 (thin solid line) in a 2048 2 system. Up 
until t = 10 all systems were evolved with the Euler update. 
In the inset is plotted the percentage difference between the 
Euler and direct updates: some error is introduced in the 
direct steps after t = 10 but at later times no increasing 
deviation from the Euler evolution is seen. 




t 

FIG. 5: Plot of scaled error per lattice site near the interface 
for a single Euler step (solid), and for a single direct step 
with ai = 3 and a 2 = (dotted with At = At 1 * 13 where from 
bottom to top A = 10~ 4 , 10" 3 , and 10~ 2 ). The scaling of 
the errors is the same as in 2] except that the errors are not 
divided by A 2 . For the two smallest A the scaling with A 2 is 
clearly seen, and so is the time independence of the error for 
the driven direct step at later times. The system size is 2048 2 
and is evolved with a Euler step with At — 0.05. 



Our single-step analysis and testing does not conclu- 
sively demonstrate that an algorithm will be reasonably 
behaved under successive steps, i.e., there is a possibility 
of accumulation of error. In Fig. we show the free en- 
ergy density for systems evolved by a direct step and com- 
pare the evolution to that obtained by the Euler method. 
It appears that the errors do not accumulate and the free 
energy decays properly as i -1 / 3 . 



C. Toward p > 2 

To go beyond the p = 2 steps with At ~ i 1 / 3 , it is 
necessary to find a stable step that satisfies Eq. i|T7|l . 
Comparing with the stability conditions, Eqs. l(2*9"jl. we 
find only marginally stable algorithms with a\ = a,2 = 
1/2 and a 3 = (3/2 - m)/(3 - to) for < m < 3. For 
to = this becomes the Crank-Nicholson method, which 
as noted before, has a fixed time step due to solvability 
considerations. However, a marginal linear iterative step 
is possible with to = 2 and a 3 — —1/2. Unfortunately, 
whether or not the marginality is a problem, the iterative 
method (given by Eyre in Q| ) fails to converge absolutely 
for these parameters. Evidently, then, it is not possible 
to construct a useful p = 3 step from the general step 
Eq. G3. 

One possible way to develop a p = 3 step is to use a 
method that is both stable and second-order accurate in 
time. For example, a two-step method that uses both 
4>t—At and cj>t to determine the updated field (f>t+At can 
be made to have no 0(At 2 ) error. A preliminary study 
of vN stability for these two-step methods indicates that 
these are a possibility. 

It is worth considering the prospect of obtaining ap-> 
oo step: according to the p classification analysis this 
would allow the natural At ~ i 2 / 3 time step. However, 
the error terms need to be strictly proportional to <9™c/> at 
each order At n . To achieve this with a one-step method 
one needs 

4>t+At - (1 - a)V 2 Mt+At =<k + aV 2 m. (40) 

Eq. I|30[) shows that this step will be linearly unstable 
when a > 1/2 (for large enough At), while for a < 1/2 
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one runs into solvability problems. At this point it seems 
unlikely that a p — > oo algorithm for the CH equation 
will be possible. 

D. Asymptotic Scaling of Field Derivatives 

In this section we derive the relations 1(33(1 and 134(1 
that provided the basis for p classification. Enough is 
known about CH dynamics that we can explicitly ana- 
lyze the leading asymptotic decay of mixed space and 
time derivatives to arbitrary order. We follow the review 
by Bray Q], and we restrict ourselves to the power-law 
scaling of these terms at sufficiently late times, where all 
observable length scales that describe the domain wall 
morphology, such as the interface curvature radii, are 
proportional to the domain size L ~ t 1 / 3 . The domain 
wall thickness w does not grow with time, so m « I 
asymptotically. However, when analyzing the fields in the 
interfacial region, defined as the locus of points within 
a distance w of a domain wall center (i.e., the surface 
<f> = 0), both length scales L and w can appear. The 
remainder of the system is referred to as the bulk. 

The scale of the chemical potential fi is proportional to 
interface curvature k due to the Gibbs-Thompson effect, 
and since k ~ 1/L 

M ~l/L~ 1/^/3. ( 41 ) 

In the bulk, the chemical potential varies smoothly and 
continuously, so a Laplacian simply brings in more pow- 
ers of L: 

V 2 M~ l/£ 3 ~ 1/t, (42) 

which implies dt4> ~ 1/t via the equation of motion JSJ. 
Now we use the relation <p — 4> eq ~ // in the bulk [l| to 
relate derivatives of <f> and A*- For example, V 2 ~ V 2 /x, 
so dt4> ~ V 2 0. Taking more time derivatives gives 

d?(f> ~ V 2 a t ""V ~ i" 2 / 3 ^" 1 ^ (43) 

Iterating this from the initial value for dt<j> gives d™(j> ~ 

i -(l/3)-(2/3)n j the firgt term in Eq ^ 

When the time derivatives act on a power of the field 
<f> 3 , the resulting expression contains the j fields and n 
time derivatives in various combinations. In this case 
the asymptotic decay comes from the single term pro- 
portional to (j>'~ 1 d™<f), which means the decay for <fp 
derivatives is the same as the j = 1 case, since the 
field <j) is order unity in the bulk. To illustrate, consider 
dftfi 3 = 6<fi(dt4>) 2 + 3(f> 2 df(j). The second term decays as 
t~ 5 / 3 as advertised, while the first term goes as (t -1 ) 2 
and is asymptotically negligible. 

Adding spatial derivatives in the bulk simply brings 
more factors of L^ 1 , so 

(V 2 ) fc 9 t ""V ~ L _a *ap _ y ~ r 2fc / 3 t-(V3)-(2/3)n 

(44) 



which gives the second term in Eq. ((3411 . 

Near interfaces, <fi changes by an amount Acj> eq in the 
amount of time, r = w/v ~ t 2 / 3 , it takes an interface to 
pass by. Therefore we get dt<fi ~ t~ 2 ^ 3 in the interfacial 
region, in contrast to dt4> ~ t . in the bulk. To de- 
termine the scaling d?(f>, consider sitting at a point just 
outside the interfacial region, in front of the moving in- 
terface. At a time 0(r) later this point will be in the 
interfacial region, so dt<p will have changed from a bulk 
to an interfacial value. This gives 

d 2 t <f> ~ (t~ 2 / 3 - t~ l )/T ~ t~ 4/3 . (45) 

Repeating this argument for higher derivatives gives 
d™4> ~ t~ 2 ™/ 3 in the interface, the first term in Eq. J33J. 

For time derivatives of cjp at the interface, we again 
get multiple terms with the various combinations of n 
time derivatives and j fields. In this case, however, every 
term contributes to the asymptotic decay. Essentially 
every time derivative, wherever it acts, brings a factor of 
t~ 2 / 3 , and these are the only factors causing the decay. 
Hence 9"<^ ~ &t<f>- Finally, adding spatial derivatives in 
the interfacial region brings factors of w _1 rather than 
L _1 , and so does not change the asymptotic decay. This 
proves the second relation in Eq. (13311 . 

IV. COMPUTATIONAL ADVANTAGE 

Having established the possibility of controlled accu- 
racy CH simulation with a growing step size At ~ t a , 
we now explore the relative computational advantage of- 
fered by such an algorithm. As described in Sec. the 
goal in such simulations is to evolve as far as close as 
possible to the scaling regime, meaning the largest pos- 
sible L(t). This means evolving until finite size effects 
enter, since stopping earlier means a smaller system size 
could be chosen. Finite size effects are expected to ap- 
pear when L(t) ~ Lot 1 / 3 is some fraction of the system 
size, so we define the simulation ending time t max by 
L(t max ) — fL S y S , or 

t max = {fL sys /L f = (fAxN/L ) 3 (46) 

where TV is the linear size of the lattice and / is a small 
constant factor. There is some arbitrariness in the defini- 
tion of the length scale L(t). We take the inverse interface 
density as our measure, that is 

using the interfacial area Ai nt from Sec. IIII Bl and its 
relation to the free energy density and surface tension 
derived therein. From our data in d = 2 we find eo ~ 
0.675, so we take L = cr/e — 1-40. 

Evolving to t max with the Euler step (or any fixed- 
size step) requires n = t max /At steps, where At is the 
step size. For our square lattice with Ax = 1 we find 
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Ato = 0.05 is close to the maximum stable value. More 
generally, one expects At ~ Ax 4 2]. Evolving to t max 
with a growing step size At ~ At" ~ dt/dn requires 

" - [ tmaX A-H- a dt ~ 1 (48) 

where a fixed-size step is used until some time to *C t max , 
and we assume to-dependent terms are negligible. 

Finally, we determine empirically the ratio of computer 
time per step j3 = Stable /TEuier- For direct steps, the FFT 
involved implies j3 ~ log N. For lattices of size 1024 2 to 
4096 2 we find ~ 2.3 ±0.1. 

Putting all this together, we find the ratio of computer 
time cost 

Euler = A(l - a)t« max = A(l - a) t fAx\ 3a r3a 
Stable (3At /3At \ L J 

(49) 

For direct steps, a — 1/3, so the relative speedup over 
Euler integration grows with the system size as N / log N. 
From At ~ Ax 4 we also see the speedup factor scaling 
as 1/Ax 3 , making stable steps an optimal choice when a 
smaller lattice spacing is desired. A p — 3 algorithm has 
a = 4/9 and offers a speedup factor of iV 4 / 3 / log N. 

We conclude by plugging in reasonably conservative 
parameter values. From Fig. 0] we see that the typical 
interfacial error for the a\ — 3, ai — direct step is about 
0.7A 2 . This is to be compared to Atfi eq — 2, the range 
in which <p varies. The choice A = 0.1 is shown in Fig. 
to give an error in the free energy density around 3% of 
the Euler value. While this seems perhaps high, we note 
that this is comparable and probably smaller than the 
error already introduced in the Euler discretization of the 
continuum CH equation due to the large lattice constant. 
It is an interesting question for future study what choice 
of Ax and A will give optimal accuracy and efficiency. 
We conclude that A — 0.1 is a reasonable choice. We 
also take a = 1/3, / = 1/10, (3 — 2.5, Ax = 1, and Lq as 
given above. These combine to give a factor 0.0387V. For 
a 1024 2 lattice the direct step is a factor 40 faster than 
the Euler method, while for a 8192 2 lattice it is a factor 
300 faster! 



V. CONCLUSIONS AND FUTURE 
DIRECTIONS 

We have seen that the general Cahn Hilliard (CH) step, 
Eq. (|10p. provides a range of linearly stable algorithms 
that prove to be gradient stable for enormous single time 
steps up to At = 10 10 . With these steps unphysical in- 
stabilities arising from the discrete implementations are 
no longer the limiting factor. Instead accuracy considera- 
tions dominate. For conserved Cahn Hilliard coarsening, 
we have analyzed and tested the accuracy scaling for sin- 
gle dynamical time steps that increase without bound 
with time as At ~ t a . We find that the errors are domi- 
nated at the order At p where they are no longer propor- 



tional to dftfi. These dominant errors restrict the growth 
of the time step to grow as At ~ i 2 (p~ 1 )/( 3 P) ; which ap- 
proaches the natural dynamical time step r <~ t 2 / 3 only as 
p — > oo. The Euler method, by contrast, is restricted to 
a constant At. This is also the case for existing implicit 
Fourier spectral algorithms. The direct steps obtained 
from Eq. I|10|l with as = 1 are linear and diagonalized 
in Fourier space, and so can be simply integrated via 
FFT's. A range of parameters, described by the shaded 
boxes in Fig. ^ are stable. These direct steps exhibit 
p = 2 and so allow At ~ t 1 / 3 , which results in speedup 
factors proportional to the linear size of the system. 

Future work in further developing these methods in- 
cludes determining possible p = 3 algorithms, for which 
At ~ t 4 / 9 is possible and the relative speedup over the 
Euler method is order N 4 / 3 / \ogN . Our preliminary 
work has shown that 0(At 2 ) accurate two-step meth- 
ods can be made unconditionally vN stable. It remains 
to test these stability predictions numerically to see if 
useful p = 3 algorithms are possible. 

It is straightforward to construct a Fourier spectral 
method integration algorithm for the stable steps ana- 
lyzed here. In fact, the numerical cost of the spectral 
method would be quite small, since the direct steps al- 
ready employ FFT's for solving the update equation. The 
primary benefit of the spectral method for unstable al- 
gorithms is that it significantly enhances the maximum 
Ato allowed by stability. It is not clear how much ben- 
efit spectral methods would bring to an already stable 
algorithm, but this should be explored. 

With the Euler step, the simulation efficiency was 
strongly dependent on Ax, leading to choosing values 
that were as large as feasible. Consequently the interface 
profile is typically poorly resolved, modifying and intro- 
ducing significant anisotropy into the surface tension. In 
contrast, the efficiency of these stable methods is much 
less dependent on the choice of lattice size, making them 
a useful tool in applications where a more accurate inter- 
face profile is desired. 

Our analysis has been for errors after a single time- 
step. If the single-step errors are small enough, the linear 
stability of bulk solutions should control the errors from 
accumulating. For the CH equation at least, our observed 
e ~ t -1 / 3 decay of the free energy, even when At ~ At 1 / 3 , 
indicates that there is no significant curvature-dependent 
modification of interfacial speeds. Nevertheless, it will 
be important to study the relationship between single- 
step errors and errors of the asymptotic scaling functions 
describing correlations to confirm this. 

We feel that our basic approach should be applicable 
in a wide variety of systems that have both nonlinear- 
ities and numerical instabilities. There are just three 
basic ingredients: i) allow for a general semi-implicit 
parametrization, following Eq. 1|10|) ; ii) check for uncon- 
ditional von Neumann (linear) stability of an individual 
update step, following Sec. Ill AI and iii) numerically test 
the vN stable algorithms for speed, accuracy, and non- 
linear stability in order to pick the best parameters for 
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further study. As long as the stability criteria are lattice 
independent, the resulting algorithms should be applica- 
ble on any regular lattice in any spatial dimension, and 
even on irregular discretizations such as used in adaptive 
mesh techniques. 
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APPENDIX A: EYRE'S THEOREM 

We repeat Eyre's stability theorem 0] here to flesh 
out the derivation for the conserved dynamics case, and 
to clarify some details of the proof. In particular, there 
are a few misleading equations in [j| that lack factors 
of the norm of the vector. More substantively, we find 
that Eyre's theorem as originally presented was slightly 
more restrictive than necessary. Note that questions of 
accuracy are not addressed in this proof, only questions 
of numerical stability. 

A central quantity in Eyre's theorem is the Hessian 
matrix 



= d 2 F 
13 dfad&j 



(Al) 



where F is the free energy and 4>i represents the field at 
the lattice site i (we consider only scalar one-component 
fields here). For free energies of interest in coarsening, 
this matrix has both positive and negative eigenvalues. 
Eyre finds a stable first-order step by splitting the free 
energy into contractive and expansive parts, F = F c + 
F E , such that F c is convex and F E is concave; that is, 
the eigenvalues of M^j, the Hessian matrix corresponding 
to F c , are strictly non- negative, and the eigenvalues of 
M E corresponding to F E are strictly non-positive, for 
any possible field configuration. 

Let X m i n < represent the lower bound for the eigen- 
values of M over all fields <p (such a bound must exist 
0), and A^ aa . < represent the upper bound on the 
eigenvalues of M E . The main result is that if 



for nonconserved dynamics or 
,SF C 



<t>t+ 



At 



AtV 1 



AtM 1 



,SF l 



4>t+&t 



(A4) 



for conserved dynamics lead to a strict non-increase of 
the free energy in time: 



(A5) 



where we have suppressed the lattice index for clarity. 
This holds unconditionally for all field configurations tfi t 
and all step sizes At > 0. Convexity of F c ensures that 
the implicit equation for 4>t+At has a unique solution. 

The energy dissipation property, along with other rea- 
sonable requirements like positivity of F, is called gra- 
dient stability by Eyre Q. While gradient stability can 
be obtained for many algorithms, such as the Euler step, 
by using a small enough At, the algorithm defined by 
Eqs. HA2|I - I|A4|I guarantees it for arbitrarily large At\ 
Even so, finding the splittings into F c and F E that lead 
to Eq. l|A2jl can be a difficult task, and the splittings, if 
they exist, may not be unique. 

Condition Eq. i|A2|l corrects the corresponding condi- 
tion in |3J , ^max — ^min ■ The current form is less restric- 
tive since X m in < 0. 

An extremely useful corollary to Eyre's theorem is that 
if the eigenvalue condition Eq. I|A2(1 is satisfied for a re- 
stricted set of fields 0, then Eq. (|A5J) still applies for all 
At provided <p t always stays within this restricted set. 
For example, <p could be field configurations with <jyf < c/)q 
for all i, for some constant 4>q. This can be useful when <fi 
is physically restricted by the dynamics, and is employed 
in the direct algorithms discussed in Sec. Ill Al 

The proof of Eq. (|A5|) relies on two inequalities 



F(cf> t+At ) - F(<j> t ) <J2< 



dF 



ix 

2 "mm 



\S4>\ 2 

(A6) 



and 



( dF E 



dF 1 



4>t+&t 



<\ta X \HV (A7) 



where 5(f>i = fcj+At ~ <f>i,t and \5<f>\ 2 = Y,i 5 4>i- These 
are simply properties of multivariable functions, and are 
derived in appendix iBl for completeness. 

Consider first nonconserved dynamics. By adding 
At[dF E /9^i]0 t+At to both sides of the equation of mo- 
tion Eq. (| A3|> one obtains 



dF 



dF 1 



OF 1 



(A8) 



\ E < l\ ■ 

"max — 2 mm 



then the field equations of motion 
SF C 



4>t+At + At- 



= 4> t - At 



5F h 



4>t+&t 



(A2) 



(A3) 



Substituting this into Eq. \A6} gives 

/ 3F E 



- 1 ix 

I 2 "Tain 



dF 1 



^)m 2 . (A!) i 
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Next use Eq. (| A7|> to complete the proof: 



F(4h+m) - F{<t> t ) < A 



: 



(A10) 

where the last inequality follows by assumption Eq. 1A2JI . 

Analyzing conserved dynamics is complicated by the 
Laplacian in the equations of motion. Consider a gen- 
eral dimensional lattice of n sites with lattice Laplacian 
(V 2 )y = Ay a symmetric n x n matrix with eigenvalues 

Ai = and A m < for all m > 1. Let uf"^ represent the 
ith component of the mth eigenvector of A, then we can 
write the Kronecker delta function as 



APPENDIX B: INEQUALITIES USED IN EYRE'S 
THEOREM 

For completeness, we re-derive Eqs. <|A6(1 and (|A7|) 
here. Consider a general function /(x) of n variables 
x = (xx,...,x n ). From the Fundamental Theorem of 
Calculus 



/(x + y)-/(x)=Y> ! ds x 

Jo 



df_ 

dxi 



(Bl) 



that is, we introduce the parameter s\ to integrate along 
the "diagonal" path from x to x + y. Similarly, we can 
write 



m— 1 



(m) (m) 



= ^2A ij A jk +u^u^ (All) 



where the pseudo-inverse A is defined by 



~ 1 
A ^ = E YI 



(m) (m) 

u\ 'u) 



(A12) 



Note that the eigenvalue Ai = corresponds to the eigen- 
vector iq — l/y/n for all i, i.e., a uniform field. Now we 
insert Eq. IjAllf) into the sum in Eq. I|A6|I and sum on k 
to get 



F(<fH+te) ~ F{4>t) < ^SfaAijAjk — 

i,j,k < P k <t>t+At 

- iA mm |£0| 2 (A13) 

where we have used ^ 5<fii = 0, which follows from the 
conservation law. Proceeding by analogy with the non- 
conserved case, we subtract AtAjk[dF E /d4>k\cpt+&t from 
both sides of the equation of motion Eq. I|A4|) to get 

dF 



At 



+ ^,A jk 

k 



0t + At 

Q F E 



dF 1 



dcj) k 



(A14) 



Substituting this into Eq. I|A13|) gives 
F(cf> t+At ) - F(&)<^W — 

i \ r 



dF 1 



dipt 



2 "min 



bjAij 



which is identical to Eq. (|A9|) except for the l/At term. 
From the definition of A and an expansion of S(f> in the 
eigenvalues u^ m > it follows that 



E< 



bjAij < 



(A16) 



so this term can be dropped from the right hand side 
of Eq. (|A15(1 and the proof follows as before to yield 
Eq. JIIJl. 



df_ 



x+siy 



dxi 



E%- 



ds^ 



d 2 f 



' dxidxj 



(B2) 



x+s 2 y 



Combining these gives the identity 

df_ 

dxi 



/(x + y)-/(x)= J> 



/ dsi / ds2 \] 
Jo Jo 



ViVj 



d 2 f 



dxidxj 



(B3) 



x+s 2 y 



Now consider the case where the eigenvalues of the matrix 
Mij = d 2 f / dxidxj are bounded from below by some 
constant X m in for all x. In this case 



E 

hi 



ViVj 



d 2 f 



dxidxj 



(B4) 



x+s 2 y 



which follows straightforwardly from an expansion of y in 
the basis of eigenvectors of M, with |y| 2 = YliVi- Thus 
we have 



/(x+y) -/(*)>£*■!£ 



dx t 



+ |A roi „|y| 2 (B5) 



where the 1/2 follows from the s integrals. Taking the 
function / to be the free energy F with x = 4>t+At and 
y = 4>t — (fit+At results in Eq. I|A6(I . 

The second inequality results from setting si = 1 in 
Eq. i|B2ll . then multiplying by yi and summing 



(A15) J2 Vi 



df_ 

dxi 



dxi 



q dxidx j 



x+sy 

(B6) 

We then use a relation similar to Eq. QB4(I . only with the 
eigenvalues of d f / dxidxj assumed to be bounded above 
by Xmax, to get 



(dl 








\dxi 


x+y 


dxi 


3 



< A„ 



(B7) 



Now we can take / = F E and x and y as before to get 
Eq. 
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